%% Aim
% Conduct statistical analysis of behavioral data 

%% Setup
clc;
clear;
tic;
root_path = pwd; 
addpath(genpath(root_path));

nrole=2;
nprice=10;
nlottery=10;

role_label={'Buyer','Seller'};
role_color_dark={[235/255 112/255 99/255],[106/255 142/255 208/255]}; 
role_color_median={[238/255 137/255 126/255],[144/255 171/255 220/255]}; 
role_color_light={[243/255 171/255 163/255],[165/255 187/255 227/255]};

%% Load data
beh_data=readtable([root_path,'\data\behavior\behavior.csv']);
sub_data=readtable([root_path,'\data\subject\subject.csv']);
nsub=height(sub_data);

%% Expected value matrix

for price=1:10
    for lottery=2:2:20
        ev.matrix.value(price,lottery/2)=nanmean(price-lottery/2);
    end
end

figure('Renderer', 'painters', 'Position', [10 10 430 350]);
    subplot(1,1,1);
    imagesc(ev.matrix.value);
    cmap_up=[linspace(1,192/256,1000)',linspace(1,0,1000)',linspace(1,0,1000)']; 
    cmap_down=[linspace(0/256,1,1000)',linspace(32/256,1,1000)',linspace(102/256,1,1000)']; 
    cmap=[cmap_down;cmap_up];
    colormap(cmap);
    set(gca,'TickDir','out');
    xlabel('Price');
    ylabel('Face Value of Lottery');
    xticks(1:1:10);
    %xticklabels({'3','4','5','6','7','8','9'});
    yticklabels({'2','4','6','8','10','12','14','16','18','20'});
    colorbar;
    %caxis([-9 9]);
    title('Relative Price');
print(1, '-dtiff', [root_path, '\figure\relative_price_matrix.tif'], '-r200');
close;

%% Probability of trading: Matrix figure (group-level)

    for irole=1:nrole
        for iprice=1:nprice
            for ilottery=1:nlottery
                bdata=beh_data(  beh_data.Role==irole & ...
                                 beh_data.Price==iprice & ...
                                 beh_data.Lottery==ilottery*2,:);
                prob_accept.matrix.group.role(irole).value(ilottery,iprice)=nanmean(bdata.Choice);
            end
        end
    end

    figure('Renderer', 'painters', 'Position', [10 10 1000 350]);
    for irole=1:nrole
        subplot(1,2,irole);
        A=prob_accept.matrix.group.role(irole).value;
        h=imagesc(A);
        cmap_up=[linspace(1,192/256,1000)',linspace(1,0,1000)',linspace(1,0,1000)']; 
        cmap_down=[linspace(0/256,1,1000)',linspace(32/256,1,1000)',linspace(102/256,1,1000)']; 
        cmap=[cmap_down;cmap_up];
        colormap(cmap);
        set(gca,'TickDir','out');
        xlabel('Price');
        ylabel('Lottery');
        yticks(1:nprice);
        xticks(1:nlottery);
        yticklabels({'2','4','6','8','10','12','14','16','18','20'});
        xticklabels({'1','2','3','4','5','6','7','8','9','10'});
        colorbar;
        caxis([0 1]);
        title(role_label{irole});
        hold on
    end
    hold off
    print(1, '-dtiff', [root_path, '\figure\behavior\choice\matrix\prob_accept_matrix_group.tif'], '-r200');
    close;        

%% Probability of trading: Descriptive results

for isub=1:nsub
    bdata=beh_data(beh_data.Subject==isub & beh_data.Role==1,:);
    sub_data.ProbBuy(isub)=nanmean(bdata.Choice);
    bdata=beh_data(beh_data.Subject==isub & beh_data.Role==2,:);
    sub_data.ProbSell(isub)=nanmean(bdata.Choice);    
end

sdata=sub_data;
A=[mean(sdata.ProbBuy);median(sdata.ProbBuy);std(sdata.ProbBuy);min(sdata.ProbBuy);max(sdata.ProbBuy);std(sdata.ProbBuy)/sqrt(nsub-1)];
B=[mean(sdata.ProbSell);median(sdata.ProbSell);std(sdata.ProbSell);min(sdata.ProbSell);max(sdata.ProbSell);std(sdata.ProbSell)/sqrt(nsub-1)];
rowNames = {'mean','median','std','min','max','se'};
colNames = {'Buy','Sell'};
prob_accept.describe = array2table([A,B],'RowNames',rowNames,'VariableNames',colNames);

% One sample t test against 50%, one tail
prob_accept.ttest{1,1}='h';
prob_accept.ttest{2,1}='p';
prob_accept.ttest{3,1}='ci';
prob_accept.ttest{4,1}='stats';
[prob_accept.ttest{1,2},prob_accept.ttest{2,2},prob_accept.ttest{3,2},prob_accept.ttest{4,2}] = ttest(sdata.ProbBuy(:),0.5);%,'Tail','left');
[prob_accept.ttest{1,3},prob_accept.ttest{2,3},prob_accept.ttest{3,3},prob_accept.ttest{4,3}] = ttest(sdata.ProbSell(:),0.5);%,'Tail','left');

disp('Probability of acceptance');
disp(prob_accept.describe);
disp('One sample ttest against 0.5, one tail');
disp(prob_accept.ttest);
disp('t stats buy');
disp(prob_accept.ttest{4,2});
disp('t stats sell');
disp(prob_accept.ttest{4,3});

%% Probability of trading: Descriptive results (Plot)

% Plot: bar
figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax=subplot(1,1,1);
ab=bar([1 2],[mean(sdata.ProbBuy) mean(sdata.ProbSell)],'EdgeColor','none','FaceColor',role_color_dark{1},'barwidth',0.6);
ab.FaceColor='flat';
ab.CData(1,:) = role_color_median{1};
ab.CData(2,:) = role_color_median{2};
hold on
er=errorbar([1],[mean(sdata.ProbBuy)],[std(sdata.ProbBuy)/sqrt(height(sdata)-1)],'Color',[0 0 0]);
er=errorbar([2],[mean(sdata.ProbSell)],[std(sdata.ProbSell)/sqrt(height(sdata)-1)],'Color',[0 0 0]);
er.LineStyle = 'none';                           
% Edit figure
xlim([0.25 2.75]);
ylim([0 1]);
set(gca,'TickDir','out');
xticklabels({'Buyer','Seller'});
ylabel('Probability of Trading');
yticks(0:0.5:1);
line([xlim],[0.5 0.5],'LineStyle','--','Color',[0 0 0]);
hold off
print(1, '-dtiff', [root_path,'\figure\behavior\choice\prob_accept_bar_group.tif'], '-r200');
close;

%% WTP and WTA: Descriptive results

beh_data.PriceLotteryEVDiff=beh_data.Price - beh_data.Lottery./2;

bdata=beh_data(beh_data.Role==1,:);
glme1 = fitglme(bdata,'Choice ~ 1 + PriceLotteryEVDiff + (1 + PriceLotteryEVDiff | Subject)','Distribution','binomial','Link','logit');
bdata=beh_data(beh_data.Role==2,:);
glme2 = fitglme(bdata,'Choice ~ 1 + PriceLotteryEVDiff + (1 + PriceLotteryEVDiff | Subject)','Distribution','binomial','Link','logit');

bdata=beh_data(beh_data.Role==1,:);
bdata.PriceLotteryEVDiff_scaled=rescale(bdata.PriceLotteryEVDiff);
glme1 = fitglme(bdata,'Choice ~ 1 + PriceLotteryEVDiff_scaled + (1 + PriceLotteryEVDiff_scaled | Subject)','Distribution','binomial','Link','logit');
bdata=beh_data(beh_data.Role==2,:);
bdata.PriceLotteryEVDiff_scaled=rescale(bdata.PriceLotteryEVDiff);
glme2 = fitglme(bdata,'Choice ~ 1 + PriceLotteryEVDiff_scaled + (1 + PriceLotteryEVDiff_scaled | Subject)','Distribution','binomial','Link','logit');

xfit=[-9:0.001:9]'; % used for identifying the indifference point

% Regression for each subject
for isub=1:nsub
    for irole=1:nrole
        % warning('Continue');
        bdata=beh_data(beh_data.Subject==isub & beh_data.Role==irole,:);
        glm = fitglm(bdata,'Choice ~ 1 + PriceLotteryEVDiff','Distribution','binomial','Link','logit');
        b(isub,1)=glm.Coefficients{1,1};
        b(isub,2)=glm.Coefficients{2,1};
        yfit{isub,irole} = glmval(b(isub,:)',xfit,'logit');
    end
end

% Identify x leading to 50% probability of trading for each subject
deltaProb=ones(nsub,1)*0.5;
for isub=1:nsub
    for j=1:18001
        if abs(yfit{isub,1}(j)-0.5)<deltaProb(isub)
            deltaProb(isub)=abs(yfit{isub,1}(j)-0.5);
            sub_data.WTP(isub)=xfit(j);            
        end
    end
end
deltaProb=ones(nsub,1)*0.5;
for isub=1:nsub
    for j=1:18001
        if abs(yfit{isub,2}(j)-0.5)<deltaProb(isub)
            deltaProb(isub)=abs(yfit{isub,2}(j)-0.5);
            sub_data.WTA(isub)=xfit(j);            
        end
    end
end
sub_data.EndowEff = sub_data.WTA - sub_data.WTP;

mean(sub_data.EndowEff)
std(sub_data.EndowEff)
[h,p,ci,stats]=ttest(sub_data.EndowEff,0,'tail','right');

% Regression for the group
for irole=1:nrole
    bdata=beh_data(beh_data.Role==irole,:);   
    glm = fitglm(bdata,'Choice ~ 1 + PriceLotteryEVDiff','Distribution','binomial','Link','logit');
    b_group(1,1)=glm.Coefficients{1,1};
    b_group(1,2)=glm.Coefficients{2,1};
    yfit_group{1,irole} = glmval(b_group(1,:)',xfit,'logit');
end

% descriptive results
sdata=sub_data; 
A=[mean(sdata.WTP);median(sdata.WTP);std(sdata.WTP);min(sdata.WTP);max(sdata.WTP)];
B=[mean(sdata.WTA);median(sdata.WTA);std(sdata.WTA);min(sdata.WTA);max(sdata.WTA)];
C=[mean(sdata.WTA-sdata.WTP);median(sdata.WTA-sdata.WTP);std(sdata.WTA-sdata.WTP);min(sdata.WTA-sdata.WTP);max(sdata.WTA-sdata.WTP)];
rowNames = {'mean','median','std','min','max'};
colNames = {'WTP','WTA','WTA_WTP'};
wtpwta.choice.describe = array2table([A,B,C],'RowNames',rowNames,'VariableNames',colNames);

% paired sample t test, one tail
wtpwta.choice.ttest{1,1}='h';
wtpwta.choice.ttest{2,1}='p';
wtpwta.choice.ttest{3,1}='ci';
wtpwta.choice.ttest{4,1}='stats';
[wtpwta.choice.ttest{1,2},wtpwta.choice.ttest{2,2},wtpwta.choice.ttest{3,2},wtpwta.choice.ttest{4,2}] = ttest(sdata.WTP(:),sdata.WTA(:));%,'Tail','left');

disp('WTP vs. WTA');
disp(wtpwta.choice.describe);
disp('paired sample ttest, one tail');
disp(wtpwta.choice.ttest);
disp('t stats');
disp(wtpwta.choice.ttest{4,2});

[h,p,ci,stats]=ttest(sdata.WTP(:),0);%,'Tail','left');
[h,p,ci,stats]=ttest(sdata.WTA(:),0);%,'Tail','left');

[h,p,ci,stats] = ttest(-sdata.WTP(:),sdata.WTA(:));%,'Tail','left');

%% WTP and WTA: Plot curves (group)

temp_role_color_light={[251/255 223/255 221/255],[229/255 235/255 247/255]};
% Plot subject curve
figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax(1)=subplot(1,1,1);
hold on;
for isub=1:nsub
    for irole=1:nrole
        plot(xfit,yfit{isub,irole}(:),'Color',temp_role_color_light{irole});
    end
end
% Plot group curve
% plot(EVxfit,mean(EVyfit)','Color',[0.4 0.2 0.6],'LineWidth',2); % based on average
for irole=1:nrole
    plot(xfit,yfit_group{1,irole},'Color',role_color_dark{irole},'LineWidth',1);
end

% Edit figure
set(ax(1), 'xlim', [-9 9]);
set(ax(1), 'ylim', [0 1]);
xticks(-9:3:9);
yticks(0:0.5:1.0);
set(ax(1),'TickDir','out');
xlabel('Price - Lottery EV');
ylabel('Probability of Trading');
% set 50% line
line(xlim,[0.5 0.5],'Color',[0 0 0],'LineStyle','--');
box on;

line([mean(sdata.WTP),mean(sdata.WTP)],[0 0.5],'Color',role_color_dark{1},'LineWidth',1,'LineStyle','--');
line([mean(sdata.WTA),mean(sdata.WTA)],[0 0.5],'Color',role_color_dark{2},'LineWidth',1,'LineStyle','--');

print(1, '-dtiff', [root_path,'\figure\behavior\choice\curve\prob_accept_curve_group.tif'], '-r200');
hold off;
close;

%%










%% Response time: Matrix figure (group-level)

    for irole=1:nrole
        for iprice=1:nprice
            for ilottery=1:nlottery
                bdata=beh_data(  beh_data.Role==irole & ...
                                 beh_data.Price==iprice & ...
                                 beh_data.Lottery==ilottery*2,:);
                rt.matrix.group.role(irole).value(ilottery,iprice)=nanmean(bdata.RT);
            end
        end
    end

    figure('Renderer', 'painters', 'Position', [10 10 1000 350]);
    for irole=1:nrole
        subplot(1,2,irole);
        A=rt.matrix.group.role(irole).value;
        h=imagesc(A);
        cmap_up=[linspace(1,192/256,1000)',linspace(1,0,1000)',linspace(1,0,1000)']; 
        cmap_down=[linspace(0/256,1,1000)',linspace(32/256,1,1000)',linspace(102/256,1,1000)']; 
        cmap=[cmap_down;cmap_up];
        colormap(cmap);
        set(gca,'TickDir','out');
        xlabel('Price');
        ylabel('Lottery');
        yticks(1:nprice);
        xticks(1:nlottery);
        yticklabels({'2','4','6','8','10','12','14','16','18','20'});
        xticklabels({'1','2','3','4','5','6','7','8','9','10'});
        colorbar;
        %caxis([1.2 3]);
        title(role_label{irole});
        hold on
    end
    hold off
    print(1, '-dtiff', [root_path, '\figure\behavior\rt\matrix\rt_matrix_group.tif'], '-r200');
    close;        

%% Response time: Descriptive results

for isub=1:nsub
    bdata=beh_data(beh_data.Subject==isub & beh_data.Role==1,:);
    sub_data.RTBuy(isub)=nanmean(bdata.RT);
    bdata=beh_data(beh_data.Subject==isub & beh_data.Role==1 & beh_data.Choice==1,:);
    sub_data.RTBuyYes(isub)=nanmean(bdata.RT);
    bdata=beh_data(beh_data.Subject==isub & beh_data.Role==1 & beh_data.Choice==0,:);
    sub_data.RTBuyNo(isub)=nanmean(bdata.RT);    
    bdata=beh_data(beh_data.Subject==isub & beh_data.Role==2,:);
    sub_data.RTSell(isub)=nanmean(bdata.RT); 
    bdata=beh_data(beh_data.Subject==isub & beh_data.Role==2 & beh_data.Choice==1,:);
    sub_data.RTSellYes(isub)=nanmean(bdata.RT);
    bdata=beh_data(beh_data.Subject==isub & beh_data.Role==2 & beh_data.Choice==0,:);
    sub_data.RTSellNo(isub)=nanmean(bdata.RT);        
end

% descriptive results

sdata=sub_data;
A=[mean(sdata.RTBuy);median(sdata.RTBuy);std(sdata.RTBuy);min(sdata.RTBuy);max(sdata.RTBuy)];
B=[mean(sdata.RTSell);median(sdata.RTSell);std(sdata.RTSell);min(sdata.RTSell);max(sdata.RTSell)];
C=[mean(sdata.RTBuyYes);median(sdata.RTBuyYes);std(sdata.RTBuyYes);min(sdata.RTBuyYes);max(sdata.RTBuyYes)];
D=[mean(sdata.RTBuyNo);median(sdata.RTBuyNo);std(sdata.RTBuyNo);min(sdata.RTBuyNo);max(sdata.RTBuyNo)];
E=[mean(sdata.RTSellYes);median(sdata.RTSellYes);std(sdata.RTSellYes);min(sdata.RTSellYes);max(sdata.RTSellYes)];
F=[mean(sdata.RTSellNo);median(sdata.RTSellNo);std(sdata.RTSellNo);min(sdata.RTSellNo);max(sdata.RTSellNo)];
G=[mean(sdata.RTBuyYes-sdata.RTBuyNo);median(sdata.RTBuyYes-sdata.RTBuyNo);std(sdata.RTBuyYes-sdata.RTBuyNo);min(sdata.RTBuyYes-sdata.RTBuyNo);max(sdata.RTBuyYes-sdata.RTBuyNo)];
H=[mean(sdata.RTSellYes-sdata.RTSellNo);median(sdata.RTSellYes-sdata.RTSellNo);std(sdata.RTSellYes-sdata.RTSellNo);min(sdata.RTSellYes-sdata.RTSellNo);max(sdata.RTSellYes-sdata.RTSellNo)];
rowNames = {'mean','median','std','min','max'};
colNames = {'Buy','Sell','BuyYes','BuyNo','SellYes','SellNo','BuyYesNo','SellYesNo'};
rt.describe = array2table([A,B,C,D,E,F,G,H],'RowNames',rowNames,'VariableNames',colNames);

% paired sample t test, one tail
rt.ttest{1,1}='h';
rt.ttest{2,1}='p';
rt.ttest{3,1}='ci';
rt.ttest{4,1}='stats';
[rt.ttest{1,2},rt.ttest{2,2},rt.ttest{3,2},rt.ttest{4,2}] = ttest(sdata.RTBuy(:),sdata.RTSell(:));%,'Tail','left');%
[rt.ttest{1,3},rt.ttest{2,3},rt.ttest{3,3},rt.ttest{4,3}] = ttest(sdata.RTBuyYes(:),sdata.RTBuyNo(:));%,'Tail','right');
[rt.ttest{1,4},rt.ttest{2,4},rt.ttest{3,4},rt.ttest{4,4}] = ttest(sdata.RTSellYes(:),sdata.RTSellNo(:));%,'Tail','right');

disp('Response time');
disp(rt.describe);
disp('Paired sample t test, one tail');
disp(rt.ttest);
disp('t stats buy vs. sell');
disp(rt.ttest{4,2});
disp('t stats buy vs. not buy');
disp(rt.ttest{4,3});
disp('t stats sell vs. not sell');
disp(rt.ttest{4,4});

mean(sdata.RTBuyYes-sdata.RTBuyNo)
std(sdata.RTBuyYes-sdata.RTBuyNo)

mean(sdata.RTSellYes-sdata.RTSellNo)
std(sdata.RTSellYes-sdata.RTSellNo)

%% Response time: Descriptive results (Plot)

figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax=subplot(1,1,1);
ab=bar([1.25 2 3 3.75],[mean(sdata.RTBuyYes) mean(sdata.RTBuyNo) mean(sdata.RTSellYes) mean(sdata.RTSellNo)],'EdgeColor','none','FaceColor',role_color_median{1},'barwidth',0.6);
ab.FaceColor='flat';
ab.CData(1,:) = role_color_median{1};
ab.CData(2,:) = role_color_light{1};
ab.CData(3,:) = role_color_median{2};
ab.CData(4,:) = role_color_light{2};
hold on
er=errorbar([1.25 2 3 3.75],[mean(sdata.RTBuyYes) mean(sdata.RTBuyNo) mean(sdata.RTSellYes) mean(sdata.RTSellNo)],[std(sdata.RTBuyYes)/sqrt(height(sdata)-1) std(sdata.RTBuyNo)/sqrt(height(sdata)-1) std(sdata.RTSellYes)/sqrt(height(sdata)-1) std(sdata.RTSellNo)/sqrt(height(sdata)-1)],'Color',[0 0 0]);
er.LineStyle = 'none';                           
% Edit figure
xlim([0.5 4.5]);
ylim([1 2.5]);
set(gca,'TickDir','out');
xticks([1.25 2 3 3.75]);
xticklabels({'Buy','Not buy','Sell','Not sell'});
ylabel('Decision Time');
yticks(1:0.5:2.5);
hold off
print(1, '-dtiff', [root_path,'\figure\behavior\rt\rt_bar_group.tif'], '-r200');
close;

%% Response time based WTP and WTA: Descriptive results

% mean response time by each value of price/lottery face value
v=unique(beh_data.PriceLotteryEVDiff);
for isub=1:nsub
    for iv=1:length(v)
        for irole=1:nrole
            bdata=beh_data(beh_data.Subject==isub & beh_data.Role==irole & beh_data.PriceLotteryEVDiff==v(iv),:);
            v_rt{irole}(isub,iv)=mean(bdata.RT);
            if height(bdata)==1
                v_rt_se{irole}(isub,iv)=0;
            else
                v_rt_se{irole}(isub,iv)=std(bdata.RT)./sqrt(height(bdata)-1);
            end
        end
    end
end
for irole=1:nrole
    mean_v_rt{irole}=mean(v_rt{irole});
    se_v_rt{irole}=std(v_rt{irole})/sqrt(height(sub_data)-1);
    median_v_rt{irole}=median(v_rt{irole});
end

% identify the v of max rt for each subject 
for irole=1:nrole
    max_v_rt{irole}=max(v_rt{irole}')';
    for isub=1:nsub
        for iv=1:length(v)
            if v_rt{irole}(isub,iv)==max_v_rt{irole}(isub)
                v_max_v_rt{irole}(isub,1)=v(iv);
            end
        end
    end
end
for irole=1:nrole
    mean_v_max_v_rt(irole)=mean(v_max_v_rt{irole});
    se_v_max_v_rt(irole)=std(v_max_v_rt{irole})/sqrt(height(sub_data));
    median_v_max_v_rt(irole)=median(v_max_v_rt{irole});
end

% descriptive results

sub_data.RTWTP(:)=v_max_v_rt{1}';
sub_data.RTWTA(:)=v_max_v_rt{2}';
sub_data.RTEndowEff = sub_data.RTWTA - sub_data.RTWTP;


sdata=sub_data; 
A=[mean(sdata.RTWTP);median(sdata.RTWTP);std(sdata.RTWTP);min(sdata.RTWTP);max(sdata.RTWTP)];
B=[mean(sdata.RTWTA);median(sdata.RTWTA);std(sdata.RTWTA);min(sdata.RTWTA);max(sdata.RTWTA)];
rowNames = {'mean','median','std','min','max'};
colNames = {'RTWTP','RTWTA'};
wtpwta.rt.describe = array2table([A,B],'RowNames',rowNames,'VariableNames',colNames);

% paired sample t test, one tail
wtpwta.rt.ttest{1,1}='h';
wtpwta.rt.ttest{2,1}='p';
wtpwta.rt.ttest{3,1}='ci';
wtpwta.rt.ttest{4,1}='stats';
[wtpwta.rt.ttest{1,2},wtpwta.rt.ttest{2,2},wtpwta.rt.ttest{3,2},wtpwta.rt.ttest{4,2}] = ttest(sdata.RTWTP(:),sdata.RTWTA(:));%,'Tail','left');

disp('RT: WTP vs. WTA');
disp(wtpwta.rt.describe);
disp('Paired sample t test, one tail');
disp(wtpwta.rt.ttest);
disp('t stats');
disp(wtpwta.rt.ttest{4,2});

[h,p,ci,stats]=ttest(sdata.RTWTP(:),0);%,'Tail','left');
[h,p,ci,stats]=ttest(sdata.RTWTA(:),0);%,'Tail','left');

%% Response time based WTP and WTA: Descriptive results (group, curve)

figure('Renderer', 'painters', 'Position', [10 10 380 350]);
ax(1)=subplot(1,1,1);
hold on;
% Plot group curve
for irole=1:nrole
    errorbar(v,mean_v_rt{irole},se_v_rt{irole},'Color',role_color_median{irole},'LineWidth',1);
end

% Edit figure
set(ax(1), 'xlim', [-10 10]);
set(ax(1), 'ylim', [1 3]);
xticks(-9:3:9);
yticks(1:1:3);
set(ax(1),'TickDir','out');
xlabel('Price - Lottery EV');
ylabel('Response Time (s)');
box on;
% Mark RT:WTP vs. WTA
for irole=1:nrole
    line([mean_v_max_v_rt(irole),mean_v_max_v_rt(irole)],ylim,'Color',role_color_median{irole},'LineWidth',1,'LineStyle','--');
end

print(1, '-dtiff', [root_path,'\figure\behavior\rt\curve\rt_curve_group.tif'], '-r200');
hold off;
close;

%%










%% Save
writetable(sub_data,[root_path,'\data\subject\subject_behavior.csv']);

%%
toc;
